C.................... XS-Strick.FOR .......................
C.. This function returns Strickland's total cross sections for a single FLIP 
C.. grid point
C.. Written by P. Richards, August 2008
      SUBROUTINE STRICKXS(IEV, !.. FLIP electron energy index
     >                  NUMPE, !.. Number of points in the FLIP energy grid
     >                 PEGRID, !.. FLIP photoelectron energy grid for 
     >                 SIGION, !.. Sum of ion state cross sections
     >                  SIGEX, !.. Sum of excited state cross sections
     >                 SIGO1D) !.. O(1D) cross section
      IMPLICIT NONE
      INTEGER I,IE,IN               !.. loop control variables
      INTEGER IDIM,EDIM             !.. Array dimension
      INTEGER NUMPE                 !.. Number of points in the FLIP energy grid
      INTEGER INEUT                 !.. Neutral ID 1=OX, 2=O2, 3=N2, 4=He
      INTEGER IEV                   !.. FLIP electron energy index
      INTEGER NUMSTATES(4)          !.. # of excited states for each species 
      PARAMETER (IDIM=99)           !.. Array dimension
      REAL THRESHE(4,IDIM)          !.. Ionization potentials
      REAL XSGRID(4,IDIM,201)       !.. cross section grid for interpolation
      REAL PEGRID(NUMPE)            !.. FLIP photoelectron energy grid for 
      REAL SIGEX(3),SIGION(3)       !.. Sum of excited and ion state cross sections
      REAL SIGI(4,201)              !.. Store total cross sections on FLIP eV grid
      REAL SIGX(4,201)              !.. Store total cross sections on FLIP eV grid
      REAL POTSUM(4),POTAV(4,201)   !.. For calculating average ionization potentials
      REAL SIGO1D                   !.. O(1D) cross section

      DATA INEUT/0/
      !.. The first time through, get the Strickland cross sections interpolated
      !.. onto the FLIP energy grid and sum the excitation cross sections
      IF(INEUT.EQ.0) THEN
         CALL GET_STRICKXS(IDIM,NUMPE,PEGRID,INEUT,NUMSTATES,
     >      THRESHE,XSGRID)
         !.. Loop to calculate the total excitation cross sections and output
         DO IN=1,INEUT
            DO IE=1,NUMPE
              SIGX(IN,IE)=0.0
              POTSUM(IN)= 0.0
              SIGI(IN,IE)=XSGRID(IN,1,IE)
              !.. Sum the excited states (NOTE: first state is ion, skip)
              DO I=2,NUMSTATES(IN)
                !.. Weighted sum of the ionization potentials 
                POTSUM(IN)=  POTSUM(IN)+THRESHE(IN,I)*XSGRID(IN,I,IE)
                SIGX(IN,IE)=SIGX(IN,IE) + XSGRID(IN,I,IE)
              ENDDO
              !.. Find the average potential energy 
              IF(SIGX(IN,IE).GT.0) POTAV(IN,IE)= POTSUM(IN)/SIGX(IN,IE)
              IF(SIGX(IN,IE).LE.0) POTAV(IN,IE)= 0.0
            ENDDO
         ENDDO
      ENDIF

      !.. Transfer saved total cross sections to FLIP cross section arrays
      DO IN=1,3
         SIGEX(IN)=SIGX(IN,IEV)
         SIGION(IN)=SIGI(IN,IEV)
      ENDDO

      SIGO1D=XSGRID(1,2,IEV)  !.. O(1D) cross section

      !..WRITE(7,'(I6,10F10.2)') IEV, PEGRID(IEV),(POTAV(I,IEV),I=1,3)

      RETURN
      END
C:::::::::::::::::::::::::::::::::::: GET_STRICKXS ::::::::::::::::::::::::::
C.. This function returns Strickland's cross sections on the FLIP photoelectron grid
C.. Written by P. Richards, August 2008
      SUBROUTINE GET_STRICKXS(IDIM,   !.. Array dimension of
     >                        EDIM,   !.. Array dimension of
     >                      PEGRID,   !.. Photoelectron energy grid 
     >                       INEUT,   !.. OUTPUT: # of neutral species
     >                   NUMSTATES,   !.. OUTPUT: # of excited states for each species
     >                     THRESHE,   !.. OUTPUT: # excitation potentials for each species
     >                      XSGRID)   !.. OUTPUT Strickland cross sections
      IMPLICIT NONE
      INTEGER I,J,K,IE              !.. loop control variables
      INTEGER IDIM,EDIM             !.. Array dimensions
      INTEGER INEUT,INEUTSAVE       !.. Neutral ID 1=OX, 2=O2, 3=N2, 4=He; saved INEUT
      INTEGER NUMPTS(4,IDIM)        !.. Number of Strickland energies
      INTEGER ISPEC                 !.. temporary neutral species indicator
      INTEGER NUMSTATES(4)          !.. # of excited states for each species 
      REAL ENERGY(4,IDIM,IDIM)      !.. Strickland energies
      REAL XSECT(4,IDIM,IDIM)       !.. Strickland cross sections
      REAL THRESHE(4,IDIM)          !.. # excitation potentials for each species 
      REAL EGRID(IDIM),XTEMP(IDIM)  !.. energy grid for interpolation
      REAL XSGRID(4,IDIM,EDIM)      !.. cross section grid for interpolation
      REAL PEGRID(EDIM)             !.. FLIP photoelectron energy grid
      REAL SCAL(4,IDIM)             !.. Strickland scaling factor; not used
      REAL EV                       !.. Temporary electron energy
      CHARACTER*65 DASHLINE         !.. for reading the dashed line
      CHARACTER*16 SPECIES          !.. temporary species marker
      CHARACTER*70 REFERENCE        !.. Reference
      CHARACTER*70 JUNK             !.. For reading any non useful line
      REAL FTXSECT                  !.. Interpolating function

      INEUTSAVE = 1 
      ISPEC=0
      REWIND(67)   !.. in case there are multiple calls
      WRITE(6,*) ' Strickland Cross sections'

      !======== Loop to read in the Strickland cross section data
 10   READ(67,'(A)',END=90) DASHLINE
         READ(67,*,END=90) INEUT         !.. species type (OX,N2,O2,He)
         
         ISPEC=ISPEC+1                   !.. count # of electronic states
         IF(INEUT.EQ.INEUTSAVE) 
     >      NUMSTATES(INEUT)=ISPEC       !.. Save # of electronic states
         IF(INEUT.NE.INEUTSAVE) ISPEC=1  !.. reset counter when species changes
         INEUTSAVE = INEUT

         READ(67,'(A)',END=90) REFERENCE   
         READ(67,'(A)',END=90) JUNK
         READ(67,*,END=90) THRESHE(INEUT,ISPEC)  !.. excitation potential
         READ(67,*,END=90) SCAL(INEUT,ISPEC)   !.. scaling factor
         READ(67,*,END=90) NUMPTS(INEUT,ISPEC)   !.. number of energies and Xsections
         READ(67,'(A)',END=90) JUNK

         !.. Now read in the energy grid and cross sections
 15      READ(67,*,END=90,ERR=15) 
     >     (ENERGY(INEUT,ISPEC,I),I=1,NUMPTS(INEUT,ISPEC))
         READ(67,*,END=90) 
     >      (XSECT(INEUT,ISPEC,I),I=1,NUMPTS(INEUT,ISPEC))

      GOTO 10
      !======== End loop to read in the Strickland cross section data

 90   CONTINUE

      !.. These loops for the interpolation onto the FLIP energy grid
      DO J=1,INEUT
        DO ISPEC=1,NUMSTATES(J)
          K=NUMPTS(J,ISPEC)
          !.. copy the Strickland values to 1-D grids for interpolation routine
          DO I=1,K
            EGRID(I)=ENERGY(J,ISPEC,I)
            XTEMP(I)=XSECT(J,ISPEC,I)
          ENDDO
          !.. Interpolation loop
          DO IE=1,EDIM
             EV=PEGRID(IE)
             !.. Interpolate vales
             XSGRID(J,ISPEC,IE)=SCAL(J,ISPEC)*
     >          FTXSECT(IDIM,1,K,EGRID,XTEMP,EV)
             IF(XSGRID(J,ISPEC,IE).LT.0) XSGRID(J,ISPEC,IE)=0
          ENDDO          
        ENDDO
      ENDDO


      RETURN
      END
C:::::::::::::::::::::::::: FTXSECT :::::::::::::::::::::::::::::::::
C.. Linear interpolation for cross sections
C.. Written by P. Richards August 2008
      REAL FUNCTION FTXSECT(IDIM,   !.. Array dimension of
     >                      EMIN,   !.. Minimum energy index
     >                      EMAX,   !.. Maximum energy index
     >                     EGRID,   !.. Energy grid
     >                     XTEMP,   !.. Xsection array
     >                        EP)   !.. Energy for desired cross section
      IMPLICIT NONE
      INTEGER IDIM,EMIN,EMAX,I1
      REAL EGRID(IDIM),XTEMP(IDIM)
      REAL EP,FACTOR

      !.. Use bisection to find the nearest energy EGRID(K) to desired
      !.. energy (EP)
      CALL BISPLTX(IDIM,EMIN,EMAX,EP,EGRID,I1)

      !.. Now interpolate to get desired production rate. Test to make
      !.. sure we haven't run off the end of the field line grid
 2    IF(I1.EQ.EMAX) THEN
         FTXSECT=0.0
      ELSE IF(XTEMP(I1).LT.0.OR.XTEMP(I1+1).LT.0) THEN
         FTXSECT=0.0
      ELSE
         FACTOR=(EP-EGRID(I1))/(EGRID(I1+1)-EGRID(I1))
         FTXSECT=XTEMP(I1)+FACTOR*(XTEMP(I1+1)-XTEMP(I1))
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::: BISPLTX ::::::::::::::::::::::::::::::::::
C--- Use bisection to find the nearest altitude ZP(K) to desired
C--- altitude (Z). Modified in Feb 93 to do both the northern and
C--- southern hemispheres.
      SUBROUTINE BISPLTX(IDIM,BEGIN,JMAX0,Z,ZP,FIRST)
      IMPLICIT REAL(A-H,L,N-Z)
      INTEGER FIRST,LAST,MID,BEGIN,IDIM
      REAL ZP(IDIM),Z

         FIRST=BEGIN
         LAST=JMAX0
         ITER=0
      !.. This section for ascending values
      IF(ZP(BEGIN).LE.ZP(JMAX0)) THEN
 10      CONTINUE
             MID=(FIRST+LAST)/2
             IF(Z.LE.ZP(MID)) THEN
                LAST=MID-1
             ELSE
                FIRST=MID+1
             ENDIF  
         IF(FIRST.LT.LAST) GO TO 10
         IF(Z.LT.ZP(FIRST)) FIRST=FIRST-1
         IF(Z.LT.ZP(BEGIN)) FIRST=BEGIN

      !.. This section for descending values
      ELSE
 110      CONTINUE
             MID=(FIRST+LAST)/2
             IF(Z.GE.ZP(MID)) THEN
                LAST=MID-1
             ELSE
                FIRST=MID+1
             ENDIF  
             ITER=ITER+1
         IF(FIRST.LT.LAST) GO TO 110
         IF(Z.GT.ZP(FIRST)) FIRST=FIRST-1
         IF(Z.LT.ZP(JMAX0)) FIRST=JMAX0
      ENDIF

 901  FORMAT(1X,4I4,1P,9E10.2)
      RETURN
      END

C:::::::::::::::::::::: GLOWSIGS :::::::::::::::::::::::::::::::::::::
C... GLOW Inelastic cross sections for electron impact
C... NOT GOOD greater than 1000 eV
      SUBROUTINE GLOWSIGS(E,SIGEX,SIGION)
      IMPLICIT NONE
      INTEGER I                    !.. loop control variable
      REAL E                       !.. electron energy
      REAL SIGEX(3),SIGION(3)      !.. Total cross section
      REAL GLOWOXX(32),GLOWOXI(32) !.. OX cross sections below 30 eV
      REAL GLOWN2X(32),GLOWN2I(32) !.. N2 cross sections below 30 eV
      REAL SIGO1D                  !.. O(1D) cross section

      !.. Cross sections at 1 eV intervals between 0.5 and 32.5 eV
      DATA GLOWOXX/2*0.0,6.250E-18,1.800E-17,2.395E-17,2.640E-17,
     > 2.690E-17,2.610E-17,2.520E-17,2.558E-17,2.918E-17,3.236E-17,
     > 3.798E-17,4.451E-17,4.890E-17,5.135E-17,5.135E-17,5.135E-17,
     > 5.020E-17,4.990E-17,4.830E-17,4.780E-17,4.610E-17,4.575E-17,
     > 4.445E-17,4.310E-17,4.270E-17,4.180E-17,4.045E-17,3.935E-17,
     > 3.87E-17,3.810E-17/
      DATA GLOWOXI/13*0.0,4.500E-19,3.800E-18,8.100E-18,1.075E-17,
     > 1.525E-17,2.442E-17,3.340E-17,4.030E-17,4.830E-17,5.510E-17,
     > 5.850E-17,6.545E-17,7.200E-17,7.500E-17,7.800E-17,8.350E-17,
     > 8.900E-17,9.20E-17,9.400E-17/
      DATA GLOWN2X/0.000E+00,0.000E+00,2.650E-16,5.300E-17,9.600E-18,
     > 2.350E-18,8.000E-19,4.885E-18,1.968E-17,4.287E-17,7.472E-17,
     > 1.069E-16,1.414E-16,1.631E-16,1.751E-16,1.786E-16,1.806E-16,
     > 1.827E-16,1.856E-16,1.875E-16,1.898E-16,1.921E-16,1.945E-16,
     > 1.955E-16,1.975E-16,2.003E-16,2.019E-16,2.046E-16,2.055E-16,
     > 2.064E-16,2.094E-16,2.100E-16/
      DATA GLOWN2I/15*0.0,6.000E-19,2.385E-18,5.335E-18,1.090E-17,
     > 1.545E-17,2.060E-17,2.700E-17,3.300E-17,3.687E-17,4.502E-17,
     > 5.370E-17,5.810E-17,6.790E-17,7.235E-17,7.680E-17,8.680E-17,
     > 9.190E-17/

      !.. Initialize cross sections
      SIGEX(1)=0.0
      SIGEX(1)=0.0
      SIGION(3)=0.0
      SIGION(3)=0.0

      !... O(lD) cross section of Henry. ).83 is GLOW factor
      SIGO1D=0.0
      IF(E.GT.1.96) SIGO1D= 0.83 * 4E-16*(1-1.96/E)**2/E

      !.. OX total excitation cross section
      IF(E.GE.30) SIGEX(1)=5.3E-16*(1.0-11.0/E)**0.5/ E**0.7
      IF(E.LT.30) SIGEX(1)=GLOWOXX(NINT(E+0.5))

      !... total excitation cross section for N2
      IF(E.GE.30) SIGEX(3)=1.95E-15*(1.0-9.0/E)**1.5/ E**0.5
      IF(E.LT.30) SIGEX(3)=GLOWN2X(NINT(E+0.5))

      !.. Get ionization cross sections
      !CALL TXSION(E,SIGION)
      !... O+ cross section 
      IF(E.GE.30) SIGION(1)=3.3E-15*(1-11.7/E)**3.2*E**(-0.6)
      IF(E.LT.30) SIGION(1)=GLOWOXI(NINT(E+0.5))

      !... N2+ cross section 
      IF(E.GE.30) SIGION(3)=8.0E-15*(1-9.0/E)**7.1*E**(-0.6)
      IF(E.LT.30) SIGION(3)=GLOWN2I(NINT(E+0.5))

      RETURN
      END
